Binding site definition on IH-AB FLAG-PURA overexpression data set

1 What was done in this report?

In this script I define PURA binding sites from the IH-AB FLAG-PURA iCLIP data set. The following steps were performed in this order:

1. Definition of Binding sites and filtering

   • PureCLIP peak calling (on merged sample1+2 sample3+4)
   • Global filter of pureclip sites (>5 cl events on pureclip site)
   • Local (genewise) filter of pureclip sites ( > 75% percentile of gene)
   • Binding site definition (5nt, Matrix approach, discard BS where center is not max pureclip, has less then 2 sites with crosslinks)
   • Reproducibility filtering (soft boundary: 0.05% percentile of cl per binding site per sample, repro in min 3 of 4 samples)
   • Comparison with old BS definition pipeline

2. Analysis of biniding behaviour

   • Analysis of bound gene types
   • Analysis of bound gene regions

2 Input

Show code

raw_path <- "/Users/melinaklostermann/Documents/projects/PURA/01_raw_data/oe-imb_koenig_2019_11_koenig_iCLIP_PURA___v2uniqMD/"

# crosslink sites in bw format
##########################################

#all 4 sample merged by Anke
bw_all_plus_path <- paste0(raw_path, "merged/bw/imb_koenig_2019_11_allsamples.v2uniqMD.duprm.plus.bw")
bw_all_minus_path <- paste0(raw_path, "merged/bw/imb_koenig_2019_11_allsamples.v2uniqMD.duprm.minus.bw")

# single samples
bw_1_plus_path <- paste0(raw_path,"bw/all/DR/imb_koenig_2019_11_sample1.v2uniqMD.duprm.plus.bw")
bw_1_minus_path <- paste0(raw_path,"bw/all/DR/imb_koenig_2019_11_sample1.v2uniqMD.duprm.minus.bw")

bw_2_plus_path <- paste0(raw_path,"bw/all/DR/imb_koenig_2019_11_sample2.v2uniqMD.duprm.plus.bw")
bw_2_minus_path <- paste0(raw_path,"bw/all/DR/imb_koenig_2019_11_sample2.v2uniqMD.duprm.minus.bw")

bw_3_plus_path <- paste0(raw_path,"bw/all/DR/imb_koenig_2019_11_sample3.v2uniqMD.duprm.plus.bw")
bw_3_minus_path <- paste0(raw_path,"bw/all/DR/imb_koenig_2019_11_sample3.v2uniqMD.duprm.minus.bw")

bw_4_plus_path <- paste0(raw_path,"bw/all/DR/imb_koenig_2019_11_sample4.v2uniqMD.duprm.plus.bw")
bw_4_minus_path <- paste0(raw_path,"bw/all/DR/imb_koenig_2019_11_sample4.v2uniqMD.duprm.minus.bw")

#
bw_all_samples_path <- "bw_4samples.RData"
bw_merges_path <- "bw_merges.RData"


# pureclip calls 
###################################
#(obtained by running pureclip on pseudo samples 1u2 and 3u4 see below)
pureclip_path <- "/Users/melinaklostermann/Documents/projects/PURA/02_R_new_pip/01-BS_def/01-BS_def_oe/PureCLIP"


# gencode annotation v31
#################################
# this is the filtered annotation as used in molitor et al
annotation <- readRDS("/Users/melinaklostermann/Documents/projects/PURA/Molitor-et-al-2022/annotation.rds")

anno_txdb <- makeTxDbFromGRanges(annotation)

output_path <- "/Users/melinaklostermann/Documents/projects/Thesis/Thesis_code/PURA/01_binding_site_definition_oePURA/"

3 Look on raw bws

sample1 is the strongest and sample2 the weakest sample. Still all samples have a very nice depth in signal.

sample	CLsites	CLevents	CLeventsPerSite
sample 1 +	4580865	10723300	2.340890
sample 1 -	4438230	10599106	2.388138
sample 2 +	3551704	7293909	2.053637
sample 2 -	3451542	7225468	2.093403
sample 3 +	3351243	7978970	2.380899
sample 3 -	3257938	7901927	2.425438
sample 4 +	4196250	9640822	2.297485
sample 4 -	4068269	9554433	2.348525
all samples +	9623431	35637001	3.703149
all samples -	9229013	35280934	3.822828

Distribution of CL sites and events in the 4 samples

4 PureCLIP peak calling

PureCLIP was performed on the bam files of 2 pseudo samples (as it can only use 2 samples at the moment): A merge of sample 1 and 2 and a merge of sample 3 and 4. The bam files of all 4 samples were obtained from Ankes pipeline. The files were merged with samtools merge. This is the bash code that was executed on the server:

Show code

# Merge
melina@gateway:~/imb_koenig_2019_11_koenig_iCLIP_PURA___v2uniqMD/bam/all/DR$ 
  conda activate

(base) melina@gateway:~/imb_koenig_2019_11_koenig_iCLIP_PURA___v2uniqMD/bam/all/DR$ 
  samtools merge sample1u2.v2uniqMD.bam imb_koenig_2019_11_sample1.v2uniqMD.duprm.bam
  imb_koenig_2019_11_sample2.v2uniqMD.duprm.bam 
  
(base)  melina@gateway:~/imb_koenig_2019_11_koenig_iCLIP_PURA___v2uniqMD/bam/all/DR$
  samtools merge sample3u4.v2uniqMD.bam imb_koenig_2019_11_sample3.v2uniqMD.duprm.bam
  imb_koenig_2019_11_sample4.v2uniqMD.duprm.bam
  
(base) melina@gateway:~/imb_koenig_2019_11_koenig_iCLIP_PURA___v2uniqMD/bam/all/DR$
  samtools index sample1u2.v2uniqMD.bam 
  
(base) melina@gateway:~/imb_koenig_2019_11_koenig_iCLIP_PURA___v2uniqMD/bam/all/DR$
  samtools index sample3u4.v2uniqMD.bam <p>

# Copy to share project
melina@gateway:/share/project/zarnack/melina/pureclip_pura_oe_pip2$ 
  cp ~/imb_koenig_2019_11_koenig_iCLIP_PURA___v2uniqMD/bam/all/DR/sample* .

# Get genome annotation
melina@gateway:/share/project/zarnack/melina/pureclip_pura_oe_pip2$
  wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_31/GRCh38.p12.genome.fa.gz
  
melina@gateway:/share/project/zarnack/melina/pureclip_pura_oe_pip2$
  gunzip GRCh38.p12.genome.fa.gz

#PureCLIP Script

#!/bin/bash <p>
#SBATCH --partition=all <p>
#SBATCH --cpus-per-task=6 <p>
#SBATCH --mem-per-cpu=32gb <p>
#SBATCH --job-name="PureCLIP_PURA_oe_pip2" <p>
#SBATCH --array=1 <p>
#SBATCH --output=pureclip_pura_o2_pip2.out    # Standard output and error log <p>
#SBATCH --error=pureclip_pura_oe_pip2.err <p>
#SBATCH --mail-user=melinaklostermann@googlemail.com <p>
#SBATCH --mail-type=ALL <p>

echo This is task \$SLURM_ARRAY_TASK_ID

conda activate
pureclip -i /share/project/zarnack/melina/pureclip_pura_oe_pip2/sample1u2.v2uniqMD.bam  
  -bai /share/project/zarnack/melina/pureclip_pura_oe_pip2/sample1u2.v2uniqMD.bam.bai  
  -i /share/project/zarnack/melina/pureclip_pura_oe_pip2/sample3u4.v2uniqMD.bam
  -bai /share/project/zarnack/melina/pureclip_pura_oe_pip2/sample3u4.v2uniqMD.bam.bai 
  -g /share/project/zarnack/melina/pureclip_pura_oe_pip2/GRCh38.p12.genome.fa 
  -o peakcalling_pura_oe_pip2_sites.bed 
  -or peakcalling_pura_oe_pip2_regions.bed 
  -nt 10


#Run job
(base) melina@gateway:~/imb_koenig_2019_11_koenig_iCLIP_PURA___v2uniqMD/PureCLIP_PURA_oe_pip2$
  sbatch jobscript_Pureclip_PURA_oe_pip2.sh <p>
(base) melina@gateway:~/imb_koenig_2019_11_koenig_iCLIP_PURA___v2uniqMD/PureCLIP_PURA_oe_pip2$
  squeue

Show code

pureclip_sites <- import(paste0(pureclip_path,"/peakcalling_pura_oe_pip2_sites.bed"), format = "bedgraph")

#clean up columns
pureclip_sites <- as.data.frame(pureclip_sites) %>% makeGRangesFromDataFrame(keep.extra.columns = T)
pureclip_sites$NA..2 <- NULL
pureclip_sites$score <- pureclip_sites$NA.
pureclip_sites$NA. <- NULL
strand(pureclip_sites) <- pureclip_sites$NA..1
pureclip_sites$NA..1 <- NULL
pureclip_sites$round_score <- round(pureclip_sites$score, digits = 1)

pureclip_sites <- keepStandardChromosomes(pureclip_sites, pruning.mode = "coarse")
pureclip_sites

GRanges object with 1895103 ranges and 2 metadata columns:
            seqnames    ranges strand |     score round_score
               <Rle> <IRanges>  <Rle> | <numeric>   <numeric>
        [1]     chr1    944315      + |   4.07434         4.1
        [2]     chr1    944323      + |   2.31949         2.3
        [3]     chr1    944344      + |   5.61261         5.6
        [4]     chr1    944348      + |   1.89390         1.9
        [5]     chr1    944354      + |   6.26169         6.3
        ...      ...       ...    ... .       ...         ...
  [1895099]     chrM       290      - |   33.9116        33.9
  [1895100]     chrM       288      - |   24.3507        24.4
  [1895101]     chrM       287      - |   62.9039        62.9
  [1895102]     chrM       286      - |   56.0864        56.1
  [1895103]     chrM       283      - |  108.5320       108.5
  -------
  seqinfo: 25 sequences from an unspecified genome; no seqlengths

5 Add annotated genes to PureCLIP sites

PureCLIP calls are loaded into R. The metadata of the resulting GRAnges object is cleaned up, as the columns get wrong names and empty columns are intruduced, when importing .bed files as GRanges objects (not shown).

Gene annotations are megred to PureCLIP sites. For this the annotation is filtered for overlapping genes. Of two overlaping genes the longer is kept.

Show code

annotation_genes <- annotation[annotation$type=="gene"] 
# sortiere anno nach laenge, damit bei uberlappenden genen das laengere (erste) behalten wird
annotation_genes <- sort(annotation_genes, by = ~ width + seqnames + start)
# remove ueberlappende gene
annotation_genes$dup <- duplicated.GenomicRanges(annotation_genes)

annotation_genes_dup_rem <- annotation_genes[annotation_genes$dup == F]

# get genes that belong to pureclip site
pureclip_sites_anno <- subsetByOverlaps(pureclip_sites, annotation_genes_dup_rem)
# if start of gene 1 = end of gene 2 keep first
idx <- findOverlaps(pureclip_sites_anno, annotation_genes_dup_rem, 
                    ignore.strand = FALSE, minoverlap = 1L, select = "first")
# remove pureclip sites that are not annotated on a gene
pureclip_sites_anno[isNA(idx)] <- NULL 
idx <- na.omit(idx)
# combine metadata of gene anno und pureclip sites
elementMetadata(pureclip_sites_anno) <- cbind(elementMetadata(pureclip_sites_anno),
                                              elementMetadata(annotation_genes_dup_rem[idx]))

# remove empty columns (empty for genes)
pureclip_sites_anno$transcript_id <- NULL
pureclip_sites_anno$transcript_type <- NULL
pureclip_sites_anno$transcript_name <- NULL
pureclip_sites_anno$transcript_support_level <- NULL
pureclip_sites_anno$havana_transcript <- NULL
pureclip_sites_anno$hgnc_id <- NULL
pureclip_sites_anno$exon_number <- NULL
pureclip_sites_anno$exon_id <- NULL
pureclip_sites_anno$protein_id <- NULL
pureclip_sites_anno$ccdsid <- NULL
pureclip_sites_anno$ont <- NULL

save(pureclip_sites_anno, 
     file=paste0(output_path,"peakcalling_pura_oe_sites_with_geneanno.RData"))
save(annotation_genes_dup_rem, 
     file=paste0(output_path,"annotation.RData"))

5.1 How many Pureclip sites are there per gene?

A problem for this normalisation are obviously genes with very few or only one pureclip call. This is the distribution of PureCLIP calls per gene:

Of all bound genes 98.45% contain more than 5 PureCLIP calls. Therefore, the calculation of the SOB on PureCLIP sites should work fine.

6 iCLIP2: loads of signal – Global and local (genewide) filter

As this data set has a very high signal depth (which is nice), also potential background signals are amplified. In order to focus on strong binding events, weaker PureCLIP sites are filtered out. A global and a local (genewide) filter are applied. The idea of the gene wide filter is that genes which were less abundant within the samples, will also have a weaker signal. A strong global filter would therefore not only filter strong binding events, but might introduce a bias towards strongly expressed genes. In the local filter one can filter more stringendly without introducing this bias. However, one has to keep in mind that the PURA homodimer in theory has 6 RNA binding domains, and its binding to a RNA might also depent on two or more adjecant binding sites.

6.1 Global filter

To choose a global filter, we can look at the distribution of the mean cl events on pureclip sites per gene and the quantiles of cl events on all pureclip sites.

As the lowest mean cl events on pureclip sites per gene are 4 and 5, a filter above 5 will already remove no background signal on these genes. As the data is from 4 replicates 5 cl events, could mean 3 samples with 1 cl and one sample with 2 cl (other distributions are also possible), which should be a minimum to trust a peak. If a filter of 5 cl events per pureclip site is applied 0 % of the pureclip sites are discarded.

6.2 Genewise filter

How many BS per gene make sense?

• Pura dimer has in theory 6 RNA Binding sites, could also bind several times the same RNA

The advantage of the genewise filter is, that it will not introduce a bias towards strongly expressed genes. However, a to stringend genewise filter might not account for more complex binding patterns. As PURA might bind at more than one site at the same time and a binding site is more than 1 nt long, the genewise filter should allow several PureCLIp sites per gene.

On the other hand for the intended following analysis in theory one BS per gene would be sufficient.

To choose a suitable percentile for filtering, percentiles in steps of 5 are calculated for all genes. These can then be plotted in a boxplot to see their distribution:

Show code

#################
#filter pureclip hits local
###############
#make purr map with different percentiles
p <- seq(from=0, to=1, by=0.05)
p_names <- purrr::map_chr(p, ~paste0(.x*100, "%"))

#write a function with purrr to calulate the different percentiles
p_funs <- purrr::map(p, ~purrr::partial(quantile, probs = .x, na.rm = TRUE)) %>% 
  purrr::set_names(nm = p_names)

#use dyplr to calculate different percentiles for each gene
quant_genewise <- pureclip_sites_anno %>% 
  as.data.frame(.) %>%
  dplyr::group_by(gene_id) %>% 
  dplyr::summarize_at("score", funs(!!!p_funs))%>% 
  as.data.frame

head(quant_genewise)

          gene_id       0%        5%      10%       15%      20%      25%
1 ENSG00000000003 0.444717 0.7685192 1.233576 1.5351030 1.809382 2.382270
2 ENSG00000000419 0.532036 0.6256254 0.697266 0.8924722 1.028882 1.400510
3 ENSG00000000971 0.570001 0.6572252 0.798329 1.1856020 1.637400 1.940510
4 ENSG00000001036 0.438706 1.0942500 1.586670 2.2371350 2.936460 3.856970
5 ENSG00000001084 0.472970 0.7554631 1.005281 1.2895080 1.944676 2.182713
6 ENSG00000001167 0.761471 0.9213473 1.317298 1.5237870 1.674410 1.920680
       30%      35%      40%      45%      50%       55%       60%       65%
1 2.955044 3.616389 4.247134 4.637863 5.418260  6.014324  6.592084  7.501156
2 1.547132 2.697808 3.326084 4.582192 5.093980  5.305564  5.821052  6.138228
3 2.096612 2.249468 2.345584 2.563632 2.705870  3.374558  4.093312  4.176944
4 4.408220 5.367835 6.505350 7.771175 9.785430 11.026100 13.308200 15.484800
5 2.491376 2.861793 3.684994 4.080322 4.528375  4.974895  5.392756  6.545365
6 2.129244 2.453282 2.532444 2.782931 3.026340  3.749743  4.165878  4.602167
        70%       75%       80%       85%       90%       95%      100%
1  8.434204  9.372795 11.839500 13.116010 17.022860 23.237190  59.02170
2  6.377186  7.351730  8.817758 11.207200 12.772940 14.816200  27.48930
3  4.350984  4.672140  5.202038  5.759970  6.458698  7.451684   9.42447
4 17.672500 20.800750 24.515000 29.131650 37.999300 53.965350 138.05600
5  7.049942  7.851968  9.007090 10.728840 14.001580 17.130815  22.56150
6  5.264186  5.953965  7.350134  8.112205  9.023546 10.631790  16.73950

Show code

# use cl event 75% quantile for filtering
quant_genewise_75 <- data.frame(gene_id = quant_genewise$gene_id, 
                                score_75_quantile = quant_genewise$`75%`)
pureclip_sites_anno<- merge(pureclip_sites_anno, 
                                          quant_genewise_75, by="gene_id")

pureclip_sites_anno <- pureclip_sites_anno[
  pureclip_sites_anno$score > pureclip_sites_anno$score_75_quantile,]


save(pureclip_sites_anno, 
     file=paste0(output_path, "peakcalling_pura_oe_sites_anno_filtered.RData"))

rm(pureclip_sites_per_gene, pureclip_sites_per_gene_plot, 
   quant_genewise, quant_genewise_75, quant_rel, quant_rel_gg)

In the end the 75% lower percentile per gene is filtered out. # Binding site definition

Binding sites are obtained, by merging all pureclip sites with a gap <= 8, and then splitting up the large BS obtained from the merging. The large binding sites are split up using a rle matrix. This approach requires afterwards to throw out Binding sites (BS) whose center is no Pureclip site or does not have the maximun Pureclip score.

Show code

###################################
# Binding site Definition
###################################
#Parameters:
#pureclip - Granges containing PureCLIP sites
# bw_plus, bw_minus - path to bw files
# windowsize - size thet Bindingsites should have (should be uneven)

Define_Binding_Sites <- function(pureclip, bw_plus, bw_minus, windowsize, out){
  
# Merge Gaps < 8 from single pureclip sites
pureclip = reduce(pureclip, min.gapwidth = 8)

#remove sites with 1 or 2 nt length
pureclip = pureclip[width(pureclip) > 2]
  
  
  bw_plus = import.bw(bw_plus, as="Rle")
  bw_minus = import.bw(bw_minus, as= "Rle")
  
  
  
  final.peaks.plus.gr <- GRanges()
  final.peaks.minus.gr <- GRanges()
  
  
  #Initialize the remaining PureCLIP CL regions to check for peaks
  remaining.regions.plus.gr <- subset(pureclip, strand == "+")
  remaining.regions.minus.gr <- subset(pureclip, strand == "-")
  
  window.radius <- (windowsize-1)/2
  while(TRUE){
    
    #No regions left to check for peaks
    if (length(remaining.regions.plus.gr) == 0 & length(remaining.regions.minus.gr) == 0){
      break
    }
    
    if (length(remaining.regions.plus.gr) != 0 ){
      #Get the raw CL counts in the remaining PureCLIP CL regions
      # returns rle list of all regions and turns it into matrix
      raw.remaining.PureCLIP.CL.regions.plus.m <- as.matrix(bw_plus[remaining.regions.plus.gr])
      
      #Identify the center of the PureCLIP CL regions (Position with max counts)
      # and store its indice
      raw.remaining.PureCLIP.CL.regions.plus.m[
        is.na(raw.remaining.PureCLIP.CL.regions.plus.m)] <- -Inf # set Na to -infinite
      max.pos.indice.plus <- max.col(raw.remaining.PureCLIP.CL.regions.plus.m, 
                                     ties.method = "first")
      
      
      #Create a peak region of xnt that is centered to the max position
      peaks.plus.gr <- remaining.regions.plus.gr
      start(peaks.plus.gr) <- start(peaks.plus.gr) + max.pos.indice.plus - 1
      end(peaks.plus.gr) <- start(peaks.plus.gr)
      peaks.plus.gr <- peaks.plus.gr + window.radius
      
      
      #Store the new peaks
      final.peaks.plus.gr <- c(final.peaks.plus.gr, peaks.plus.gr)
      
      #Remove the peaks from the CL regions to search for additional peaks
      #Excise additionally x nucleotides up and downstream
      peaks.plus.grl <- as(peaks.plus.gr+window.radius, "GRangesList")
      
      remaining.regions.plus.gr <- unlist(psetdiff(remaining.regions.plus.gr, peaks.plus.grl))
    }
    if (length(remaining.regions.minus.gr) != 0 ){
      #Get the raw CL counts in the remaining PureCLIP CL regions
      # returns rle list of all regions and turns it into matrix
      raw.remaining.PureCLIP.CL.regions.minus.m <- as.matrix(
        bw_minus[remaining.regions.minus.gr])
      
      #Identify the center of the PureCLIP CL regions (Position with max counts) 
      #and store its indice
      raw.remaining.PureCLIP.CL.regions.minus.m[
        is.na(raw.remaining.PureCLIP.CL.regions.minus.m)] <- -Inf
      max.pos.indice.minus <- max.col(raw.remaining.PureCLIP.CL.regions.minus.m, ties.method = "last")
      
      #Create a peak region of xnt that is centered to the max position
      peaks.minus.gr <- remaining.regions.minus.gr
      start(peaks.minus.gr) <- start(peaks.minus.gr) + max.pos.indice.minus - 1
      end(peaks.minus.gr) <- start(peaks.minus.gr)
      peaks.minus.gr <- peaks.minus.gr + window.radius
      
      #Store the new peaks
      final.peaks.minus.gr <- c(final.peaks.minus.gr, peaks.minus.gr)
      
      #Remove the peaks from the CL regions to search for additional peaks
      #Excise additionally x nucleotides up and downstream
      peaks.minus.grl <- as(peaks.minus.gr+window.radius, "GRangesList")
      
      remaining.regions.minus.gr <- unlist(psetdiff(remaining.regions.minus.gr,
                                                    peaks.minus.grl))
    }
  }
  export(final.peaks.plus.gr, 
         con= paste(out,"_plus.bed", sep = ""), 
         format = "bed")
  export(final.peaks.minus.gr, 
         con= paste(out,"_minus.bed", sep = ""), 
         format = "bed")
  save(final.peaks.minus.gr, file= paste(out,"_plus.RData", sep = ""))
  save(final.peaks.plus.gr, file= paste(out,"_minus.RData", sep = ""))
  
  returnlist <- list(peaks.minus = final.peaks.minus.gr, peaks.plus = final.peaks.plus.gr)
  return(returnlist)
}

pureclip_sites_anno <- makeGRangesFromDataFrame(pureclip_sites_anno,
                                                             keep.extra.columns = TRUE)

binding_sites <- Define_Binding_Sites(pureclip = pureclip_sites_anno, 
                                      bw_plus = bw_all_plus_path, bw_minus = bw_all_minus_path,
                                      windowsize = 5, out = paste0(output_path,"Binding_site_windows_5nt" ))

6.3 Throw out BS where not pureclip site or not max pureclip site

Show code

load(paste0(output_path,"Binding_site_windows_5nt_plus.RData"))
load(paste0(output_path,"Binding_site_windows_5nt_minus.RData"))

############################
# Keep only BS with PureCLIP center
############################
# get all BS
binding_sites <- c(final.peaks.minus.gr, final.peaks.plus.gr)
#get centers
BS_centers <- binding_sites - 2

#keep only overlaps with pureclip sites
pureclip_sites_anno <- makeGRangesFromDataFrame(pureclip_sites_anno, 
                                                             keep.extra.columns = TRUE)
binding_sites_center_PS <- binding_sites[queryHits(findOverlaps(
  BS_centers, pureclip_sites_anno))]



###########################
# Keep only BS with max PureCLIP site at center
##########################
# get bw rles
bw_plus_rle <- import.bw(bw_all_plus_path, as="Rle")
bw_minus_rle <- import.bw(bw_all_minus_path, as="Rle")

# split BS by strand
binding_sites_center_PS_plus <- binding_sites_center_PS[strand(binding_sites_center_PS)=="+"]
binding_sites_center_PS_minus <- binding_sites_center_PS[strand(binding_sites_center_PS)=="-"]

# make matrix of BS
binding_sites_center_PS_plus_m <- as.matrix(bw_plus_rle[binding_sites_center_PS_plus])
binding_sites_center_PS_minus_m <- as.matrix(bw_minus_rle[binding_sites_center_PS_minus])

# calc max for each BS (one BS is one row in the matrix)
max_BS_plus <- apply(binding_sites_center_PS_plus_m,1,max)
max_BS_minus <- apply(binding_sites_center_PS_minus_m,1,max)

# subset for center = max
binding_sites_center_PSmax_plus <- binding_sites_center_PS_plus[
  max_BS_plus == binding_sites_center_PS_plus_m[,3]] 
binding_sites_center_PSmax_minus <- binding_sites_center_PS_minus[
  max_BS_minus == binding_sites_center_PS_minus_m[,3]]



###########################
# Keep only BS with at least 2  crosslink sites
############################
binding_sites_center_PSmax_plus_m <- as.matrix(bw_plus_rle[binding_sites_center_PSmax_plus])
binding_sites_center_PSmax_minus_m <- as.matrix(bw_minus_rle[binding_sites_center_PSmax_minus])

crosslink_sites_plus <- apply(binding_sites_center_PSmax_plus_m, 1, function(x) 5-sum(x == 0)) 
crosslink_sites_minus <- apply(binding_sites_center_PSmax_minus_m, 1, function(x) 5-sum(x == 0)) 


binding_sites_center_PSmax_plus_2cl <- binding_sites_center_PSmax_plus[crosslink_sites_plus > 1]
binding_sites_center_PSmax_minus_2cl <- binding_sites_center_PSmax_minus[crosslink_sites_minus > 1]


###########################
# BS definition steps table
##########################
BS_steps_list <- list(final.peaks.plus.gr, final.peaks.minus.gr, 
                      binding_sites_center_PS_plus, binding_sites_center_PS_minus, 
                      binding_sites_center_PSmax_plus, binding_sites_center_PSmax_minus,
                      binding_sites_center_PSmax_plus_2cl, binding_sites_center_PSmax_minus_2cl )
n_BS_steps <- sapply(BS_steps_list, NROW)
names_BS_steps <- c("definded BS plus", "defined BS minus", 
                    "PureCLIP center plus", "PureClip center minus",
                    "PureCLIP max center plus", "PureCLIP max center minus",
                    "with 2 CL sites plus", "with 2CL sites minus")

table_BS_steps <- cbind(names_BS_steps, n_BS_steps)
kable(table_BS_steps)
save(table_BS_steps, file = paste0(pureclip_path,"BS_def_steps_oe.RData"))

##########################
# final Binding sites
#########################
binding_sites_final <- c(binding_sites_center_PSmax_plus_2cl, 
                         binding_sites_center_PSmax_minus_2cl)%>%sort

# add center sob, pureclip score and cl events
# get pureclip sites that overlap with center peaks
pureclip_sites_overlapping_BS_center <- subsetByOverlaps(pureclip_sites_anno,
                                                         (binding_sites_final-2)) %>% sort
# when both Granges are sorted, they have the same order, 
# metadata can be moved columnwise without matching
# binding_sites_final$gene_id <- pureclip_sites_overlapping_BS_center$gene_id
binding_sites_final$SOB <- pureclip_sites_overlapping_BS_center$SOB
binding_sites_final$score <- pureclip_sites_overlapping_BS_center$score
binding_sites_final$cl_events_center <- pureclip_sites_overlapping_BS_center$cl_events

save(binding_sites_final, file=paste0(pureclip_path,"Binding_site_windows_5nt_final.RData"))

7 Sample Reproducibility and Reproducibility filtering

Binding sites are demed reprodrucible, if the they are found in at least three of the four binding sites. The weakest BS per sample with the 0.05 percentile of CL events per BS are not considered for the reproducibility.

per0.05	per0.1	per0.2	variable
5	7	9	clp_rep1
3	4	6	clp_rep2
3	5	7	clp_rep3
5	6	8	clp_rep4

cutoff for the four samples, at 5%, 10% and 20% quantile

Show code

#########################
# mark samplewise the lowest 20% percentile of cl per BS 
########################
binding_sites_cl_samples$cl1_overThresh <- 
  binding_sites_cl_samples$clp_rep1 > quantiles_samples[1,3]

binding_sites_cl_samples$cl2_overThresh <- 
  binding_sites_cl_samples$clp_rep2 > quantiles_samples[2,3]

binding_sites_cl_samples$cl3_overThresh <- 
  binding_sites_cl_samples$clp_rep3 > quantiles_samples[3,3]

binding_sites_cl_samples$cl4_overThresh <- 
  binding_sites_cl_samples$clp_rep4 > quantiles_samples[4,3]

# in how many samples is the number of cls per bs over the threshold?
binding_sites_cl_samples$nsamples_overThresh <-
  rowSums(cbind(binding_sites_cl_samples$cl1_overThresh,
                binding_sites_cl_samples$cl2_overThresh,
                binding_sites_cl_samples$cl3_overThresh,
                binding_sites_cl_samples$cl4_overThresh))

# keep only BS with cls over threshold in at least 3 of the 4 samples
binding_sites_repro <- 
  binding_sites_cl_samples[binding_sites_cl_samples$nsamples_overThresh>=3,] %>%
  makeGRangesFromDataFrame(., keep.extra.columns = T)

#save(binding_sites_repro, file=paste0(output_path,"final_repro_BS_oe.RData" ))

As can be seen in the reproducibility matrix the BS of all samples correlate well together. The ECDF of the coefficent variation also shows a good reporducibility as a ECDF of 1 is reached at around 0.7 (1 is proposed as a cutoff). The cv depends on the standard deviation. Probably it is more meaningfull to look at the standard error of the experiment> As SE = SD /(number of observations per sample)^-0.5 for 4 samples the SE is half of the SD.

8 Bound Genes and regions

8.1 assign bound gene type

• many binding sites overlap with multiple genes (non codings!)
• in IGV it look like rather the non-coding RNA than the protein-coding RNA is bound in that case
• we therefore use a hieracry on gene types: lncRNA > miRNA > miscRNA > snRNA > snoRNA > protein-coding

Show code

############################
# annote binding sites with the bound ( = overlapping) gene
################################
#  overlaps are resolved by:
# 1. for two genes of different type by the type hieracy
# 2. for to protein coding genes: in this chunk one protein coding gene is assigned at random in the next chunk the region hieracry is used to make this decision and the gene info is changed then

binding_sites <- binding_sites_repro 

annotation_gene <- annotation[annotation$type=="gene"]

# overlaps of bs with any gene
ol <- findOverlaps(annotation_gene, binding_sites, ignore.strand=F)

# index binding sites
binding_sites$idx <-  1:NROW(binding_sites)

# get all possible regions overlapping with crosslinked nucleotide in a temporary file
binding_sites_temp1 <- binding_sites[subjectHits(ol)]
binding_sites_temp1 <- sortSeqlevels(binding_sites_temp1)
elementMetadata(binding_sites_temp1) <-  c(elementMetadata(binding_sites_temp1), elementMetadata(annotation_gene[queryHits(ol), c("gene_id", "gene_type", "gene_name")]))

table(binding_sites_temp1$gene_type)

        lncRNA          miRNA       misc_RNA protein_coding       ribozyme 
          1944             27             46         113430             11 
        scaRNA         snoRNA          snRNA 
             3             14             48 

Show code

# chose for each binding site a gene from the highest possible hierarchy
# if two genes from the same type would be present the first one is chosen (randomly)
binding_sites_temp2 <- binding_sites_temp1 %>%  
  as.data.frame(.) %>% # this file contains several entries for 1 binding site overlapping with more than one gene
  group_by(idx) %>% # the idx colum is an index of the binding sites
  arrange(factor(gene_type, levels = c("lncRNA", "miRNA", "miscRNA", "snRNA", "snoRNA", "protein_coding" )), .by_group = T) %>% # arrange by hierarcy
  dplyr::slice(1) # choose randomly  the first gene (the chosen gene is changed in the next chunk for protein coding genes)

binding_sites_temp2 <- makeGRangesFromDataFrame(binding_sites_temp2, keep.extra.columns = T)

table(binding_sites_temp2$gene_type)

        lncRNA          miRNA       misc_RNA protein_coding         snoRNA 
          1933             27             14         106845             13 
         snRNA 
            36 

Show code

# non coding rnas "fake" region non coding
bs_temp_non_cod <- binding_sites_temp2[binding_sites_temp2$gene_type!= "protein_coding",] %>% makeGRangesFromDataFrame(., keep.extra.columns = T)
bs_temp_non_cod$type <- "non_coding"

8.2 Binding site region by hiercary approach

• Overlap Bs with mutiple regions and choose one with highest priority
• hieracry: “three_prime_UTR” > “five_prime_UTR” > “CDS” > “intron”
• only done for protein-coding genes

Show code

# classify only protein coding genes, regions only from bound genes
anno_regions <- annotation[annotation$gene_type == "protein_coding" & annotation$type %in% c("three_prime_UTR", "five_prime_UTR", "CDS") ]

#get introns
introns <- intronsByTranscript(anno_txdb, use.names=T) %>% unlist(.)
introns$type = "intron"

anno_regions <- c(anno_regions, introns)

#overlaps of binding sites with protein coding genes
binding_sites_temp_pc <- binding_sites_temp2[binding_sites_temp2$gene_type=="protein_coding"]


ol_bs <- findOverlaps(anno_regions, binding_sites_temp_pc-2, ignore.strand = FALSE)


# get all possible regions overlapping with bs
binding_sites_temp_pc2 <- binding_sites_temp_pc[subjectHits(ol_bs)]
binding_sites_temp_pc2 <- sortSeqlevels(binding_sites_temp_pc2)

# get possible region and transcript annotation of bs
anno_bs <- elementMetadata(anno_regions[queryHits(ol_bs)]) %>%
  as.data.frame(.) %>%
  dplyr::select(c("type", "gene_id"))

# and add to temp binding sites
colnames(anno_bs) <- c("type", "gene_id.from_region")
elementMetadata(binding_sites_temp_pc2) <- c(elementMetadata(binding_sites_temp_pc[subjectHits(ol_bs)]), anno_bs)


# choose type highest in hiearcy
binding_sites_temp_pc3 <- as.data.frame(binding_sites_temp_pc2) %>%
    mutate(., gene_id = gene_id.from_region ) %>% # this throws out regions that overlap with non-coding RNAs, with this line no binding sites in introns!!!
  group_by(idx) %>% # specific id per binding site
  arrange(factor(type, levels = c("three_prime_UTR", "five_prime_UTR", "CDS", "intron")), .by_group = T) %>%
  dplyr::slice(1) # arrange by region, select first, transcript info is ignored here

binding_sites_temp_pc3 <- makeGRangesFromDataFrame(binding_sites_temp_pc3, keep.extra.columns = T)
sortSeqlevels(binding_sites_temp_pc3)

GRanges object with 103555 ranges and 16 metadata columns:
           seqnames              ranges strand |     score  clp_rep1  clp_rep2
              <Rle>           <IRanges>  <Rle> | <numeric> <numeric> <numeric>
       [1]     chr1     1022219-1022223      + |   16.6259        14         7
       [2]     chr1     1022234-1022238      + |   16.1704         9        10
       [3]     chr1     1022243-1022247      + |   12.6694         7         8
       [4]     chr1     1478169-1478173      + |   34.1607        26        20
       [5]     chr1     1478176-1478180      + |   36.6225        11        20
       ...      ...                 ...    ... .       ...       ...       ...
  [103551]     chrX 154653297-154653301      - |   57.4722        39        25
  [103552]     chrX 154653344-154653348      - |   60.7973        43        37
  [103553]     chrX 154653366-154653370      - |   39.7690        42        28
  [103554]     chrX 154653375-154653379      - |   40.9502        18        25
  [103555]     chrX 155061797-155061801      - |   15.4820        12         7
            clp_rep3  clp_rep4 cl1_overThresh cl2_overThresh cl3_overThresh
           <numeric> <numeric>      <logical>      <logical>      <logical>
       [1]         9        10           TRUE           TRUE           TRUE
       [2]        12        12          FALSE           TRUE           TRUE
       [3]         9        11          FALSE           TRUE           TRUE
       [4]        15        18           TRUE           TRUE           TRUE
       [5]        22        14           TRUE           TRUE           TRUE
       ...       ...       ...            ...            ...            ...
  [103551]        18        32           TRUE           TRUE           TRUE
  [103552]        31        37           TRUE           TRUE           TRUE
  [103553]        35        29           TRUE           TRUE           TRUE
  [103554]        15        20           TRUE           TRUE           TRUE
  [103555]        10         8           TRUE           TRUE           TRUE
           cl4_overThresh nsamples_overThresh       idx         gene_id
                <logical>           <numeric> <integer>     <character>
       [1]           TRUE                   4         1 ENSG00000188157
       [2]           TRUE                   3         2 ENSG00000188157
       [3]           TRUE                   3         3 ENSG00000188157
       [4]           TRUE                   4         6            <NA>
       [5]           TRUE                   4         7            <NA>
       ...            ...                 ...       ...             ...
  [103551]           TRUE                   4    108863 ENSG00000126890
  [103552]           TRUE                   4    108864 ENSG00000126890
  [103553]           TRUE                   4    108865 ENSG00000126890
  [103554]           TRUE                   4    108866 ENSG00000126890
  [103555]          FALSE                   3    108867 ENSG00000182712
                gene_type   gene_name            type gene_id.from_region
              <character> <character>        <factor>         <character>
       [1] protein_coding        AGRN          CDS        ENSG00000188157
       [2] protein_coding        AGRN          CDS        ENSG00000188157
       [3] protein_coding        AGRN          CDS        ENSG00000188157
       [4] protein_coding      ATAD3B          intron                <NA>
       [5] protein_coding      ATAD3B          intron                <NA>
       ...            ...         ...             ...                 ...
  [103551] protein_coding       CTAG2 CDS                 ENSG00000126890
  [103552] protein_coding       CTAG2 CDS                 ENSG00000126890
  [103553] protein_coding       CTAG2 CDS                 ENSG00000126890
  [103554] protein_coding       CTAG2 CDS                 ENSG00000126890
  [103555] protein_coding        CMC4 three_prime_UTR     ENSG00000182712
  -------
  seqinfo: 25 sequences from an unspecified genome; no seqlengths

Show code

# merge BS from non-coding RNAs and binding sites with specific regions
binding_sites_with_regions <- c(binding_sites_temp_pc3, bs_temp_non_cod)


# some protein coding genes not have a region, because no trustworth transcript is annotated in this position of the gene
bs_without_region <- binding_sites[-queryHits(findOverlaps(binding_sites, binding_sites_with_regions))]

9 Output

Show code

saveRDS(binding_sites_with_regions,  paste0(output_path , "bs_oe.rds"))

10 Session Info

Show code

sessionInfo()

R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur ... 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] Gviz_1.42.1            ggpubr_0.6.0           dplyr_1.1.3           
 [4] GenomicFeatures_1.50.4 AnnotationDbi_1.60.2   Biobase_2.58.0        
 [7] rtracklayer_1.58.0     GenomicRanges_1.50.2   GenomeInfoDb_1.34.9   
[10] IRanges_2.32.0         S4Vectors_0.36.2       BiocGenerics_0.44.0   
[13] ggplot2_3.4.4          knitr_1.45            

loaded via a namespace (and not attached):
  [1] colorspace_2.1-0            deldir_1.0-9               
  [3] ggsignif_0.6.4              rjson_0.2.21               
  [5] biovizBase_1.46.0           htmlTable_2.4.2            
  [7] XVector_0.38.0              base64enc_0.1-3            
  [9] dichromat_2.0-0.1           rstudioapi_0.15.0          
 [11] farver_2.1.1                bit64_4.0.5                
 [13] fansi_1.0.5                 xml2_1.3.5                 
 [15] codetools_0.2-19            cachem_1.0.8               
 [17] polyclip_1.10-6             Formula_1.2-5              
 [19] jsonlite_1.8.7              Rsamtools_2.14.0           
 [21] broom_1.0.5                 cluster_2.1.4              
 [23] dbplyr_2.4.0                png_0.1-8                  
 [25] ggforce_0.4.1               compiler_4.2.2             
 [27] httr_1.4.7                  backports_1.4.1            
 [29] Matrix_1.5-4.1              fastmap_1.1.1              
 [31] lazyeval_0.2.2              cli_3.6.1                  
 [33] tweenr_2.0.2                htmltools_0.5.7            
 [35] prettyunits_1.2.0           tools_4.2.2                
 [37] gtable_0.3.4                glue_1.6.2                 
 [39] GenomeInfoDbData_1.2.9      reshape2_1.4.4             
 [41] rappdirs_0.3.3              Rcpp_1.0.11                
 [43] carData_3.0-5               vctrs_0.6.4                
 [45] Biostrings_2.66.0           xfun_0.41                  
 [47] stringr_1.5.0               lifecycle_1.0.3            
 [49] ensembldb_2.22.0            restfulr_0.0.15            
 [51] rstatix_0.7.2               XML_3.99-0.15              
 [53] MASS_7.3-60                 zlibbioc_1.44.0            
 [55] scales_1.2.1                BSgenome_1.66.3            
 [57] VariantAnnotation_1.44.1    ProtGenerics_1.30.0        
 [59] hms_1.1.3                   MatrixGenerics_1.10.0      
 [61] parallel_4.2.2              SummarizedExperiment_1.28.0
 [63] AnnotationFilter_1.22.0     RColorBrewer_1.1-3         
 [65] yaml_2.3.7                  curl_5.1.0                 
 [67] memoise_2.0.1               gridExtra_2.3              
 [69] UpSetR_1.4.0                biomaRt_2.54.1             
 [71] rpart_4.1.21                latticeExtra_0.6-30        
 [73] stringi_1.7.12              RSQLite_2.3.2              
 [75] BiocIO_1.8.0                plotrix_3.8-3              
 [77] checkmate_2.3.0             filelock_1.0.2             
 [79] BiocParallel_1.32.6         rlang_1.1.2                
 [81] pkgconfig_2.0.3             matrixStats_1.0.0          
 [83] bitops_1.0-7                evaluate_0.23              
 [85] lattice_0.22-5              purrr_1.0.2                
 [87] labeling_0.4.3              GenomicAlignments_1.34.1   
 [89] htmlwidgets_1.6.2           bit_4.0.5                  
 [91] tidyselect_1.2.0            plyr_1.8.9                 
 [93] magrittr_2.0.3              R6_2.5.1                   
 [95] generics_0.1.3              Hmisc_5.1-1                
 [97] DelayedArray_0.24.0         DBI_1.1.3                  
 [99] pillar_1.9.0                foreign_0.8-85             
[101] withr_2.5.2                 KEGGREST_1.38.0            
[103] abind_1.4-5                 RCurl_1.98-1.13            
[105] nnet_7.3-19                 tibble_3.2.1               
[107] crayon_1.5.2                car_3.1-2                  
[109] interp_1.1-4                utf8_1.2.4                 
[111] BiocFileCache_2.6.1         rmarkdown_2.25             
[113] jpeg_0.1-10                 progress_1.2.2             
[115] data.table_1.14.8           blob_1.2.4                 
[117] digest_0.6.33               tidyr_1.3.0                
[119] munsell_0.5.0